Flow analysis from cumulants: a practical guide 
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We have recently proposed a new method of flow analysis, based on a cumulant expansion of 
multiparticle azimuthal correlations. Here, we describe the practical implementation of the method. 
The major improvement over traditional methods is that the cumulant expansion eliminates order 
by order correlations not due to flow, which are often large but usually neglected. 



... I. INTRODUCTION 

' The measurement of the azimuthal distributions of outgoing particles with respect to the reaction plane in noncentral 

^vq ' heavy ion collisions — the flow analysis — is an important probe of the interaction region of the collision [Q . In particular, 

I . it has raised much interest at ultrarelativistic energies Q where it may signal the formation of a quark-gluon plasma 

O In addition, combining flow and two-particle interferometry results yields a three-dimensional picture of the 

. emitting source [^,^. Therefore, accurate flow measurements are highly needed. 

' Azimuthal distributions are characterized by the Fourier coefficients 
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m(0-*H)\ = (cosn(0-<ffl)), (1) 
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^ ■ where </> is the azimuthal angle of an emitted particle in the laboratory frame, is the azimuth of the reaction plane, 
' and angular brackets denote a statistical average over many particles and events. Ideally, w„ should be measured for 
various particles as a function of their transverse momentum px and rapidity y ("differential" flow). The first two 
harmonics vi and V2 are the so-called directed and elliptic flows, respectively. 

Since the azimuth in a given event is unknown, the coefficients u„ are extracted from the azimuthal correlations 
between outgoing particles. The underlying idea is that the correlation of every particle with the reaction plane 
induces correlations between the particles. Standard methods extract flow from two-particle azimuthal correlations, 

■ either directly [Q, or through the correlation between two "subevents" [^,^. However, the correlation between 
two given particles is not only due to flow, and the other sources of correlation — as, e.g., quantum Bose-Einstein 

I , . effects, momentum conservation, resonance decays, jets — may dominate the measured signal, especially for peripheral 
' collisions, thereby spoiling the validity of the analysis The impact of "nonflow" correlations on the flow 

^ , analysis might tentatively be minimized: cuts in phase space can be used to avoid the influence of quantum effects and 
' resonance decays, while the contribution of momentum conservation to the measured correlation can be calculated and 
^ , subtracted |]ri|,^. However, the various recipes require some a priori knowledge of nonflow correlations; furthermore 
it is necessary to assume that all sources of such correlations are known and accounted for, which may not be true. 
To remedy the contamination from nonflow correlations, which amounts to systematic uncertainties on the flow 

■ values, we have introduced new methods of flow analysis, based on a cumulant expansion of multiparticle azimuthal 
■ " ' correlations jljjlj]. The principle of the methods is that when cumulants of higher order are considered, the relative 

contribution of nonflow effects, and thus the corresponding systematic error, decreases. Recently, this cumulant 
expansion has been successfully applied by the STAR Collaboration at RHIC (| jl^ . 

More precisely, the cumulant of 2fc-particle azimuthal correlations, which we denote c„{2fc} (where n is the Fourier 
harmonic and 2k is an even integer, in practice 2, 4 or 6), is a quantity built with all the measured azimuthal 
correlations up to order 2k, i.e., the {exp[in{(f>i + (f>2 + ■ ■ ■ + (f>k' — (pk'+i — ■ ■ ■ — (j)k'+k")]), with k' + k" < 2k. The 
key feature of the cumulant is that it eliminates the contribution of lower order correlations, so that only the genuine 
2fc-particle correlation remains. Flow, which is essentially a collective effect, gives a contribution to the cumulant 
proportional to vf^^. The remaining contribution, from 2fc-particle nonflow correlations, scales as N^~'^^, where N is 
the total multiplicity of particles emitted in an event [^3|. Therefore, the flow dominates if 

Stated differently, if the cumulant c„{2fc} is much larger than N^^^"^^ , it is dominated by flow, and therefore yields 
an estimate of w„, which we denote i;„{2fc}, with a systematic error (due to unknown nonflow correlations) of order 
0{{NvnY~'^^)- When k increases, the systematic error decreases: this is why the estimate w„{4} derived from the 
fourth cumulant is a priori more accurate than Vn{2}, which is in fact the value given by two-particle methods. 
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In the following, we describe the practical implementation of the method, referring the reader to Refs. [g^ 14 
for further theoretical justifications. The first step consists in deriving a global measurement of t;„, integrated over 
some phase-space region, typically a detector acceptance, for a given centrality class (Sec. |l|): this is equivalent 
to reconstructing the reaction plane and obtaining the "event plane resolution" in the subevcnt method g|. This 
"integrated" flow serves as reference for the differential flow analysis discussed in Sec. III. An important feature of 



our method is that it automatically takes into account azimuthal inhomogeneities in the acceptance of the detector. 



In the case of a detector with only partial azimuthal coverage, minor modifications occur, which are given in Sec. IV 



II. INTEGRATED FLOW 

Consider a data set of iVovts events of approximately the same centrality, recorded in a run with a constant detector 
acceptance. In this section, we explain how estimates of the flow, integrated over the detector acceptance, can be 
obtained from cumulants of 2-, 4-, and 6-particle correlations. 

We denote by (f>j the azimuths of the outgoing particles with respect to a fixed direction in the laboratory. The 
various quantities of interest are constructed from the real- valued generating function [p^ 

M M 

Gn{z) = n [l + ^ i'*''"^' + ze--*0] = n [l + ^ {2xcos{n^,) + 2ysin(n0,))] , (3) 

where the product runs over M particles detected in a single event and z = x + iy is an arbitrary complex number. 
This generating function has no physical meaning in itself, but after averaging over events, the coefficients of its 
expansion in powers of z and z* = x — iy yield multiparticle azimuthal correlations of arbitrary orders. In practice, 
the variable z corresponds to interpolation points used to estimate the various quantities encountered in the analysis, 
as will be explained shortly. 

As in the standard flow analysis, a weight Wj is attributed to particle j, which is a function of particle type, transverse 
momentum, and rapidity. Naturally, the integrated flow obtained from this generating function will be weighted by 
w, I.e. weight must be chosen so as to maximize the effects of flow relative to statistical 

fluctuations. As we shall see below, this is achieved by maximizing the dimensionless quantity Xn = Vny^ M/ {w"^) 
(this quantity also characterizes the event plane resolution in the standard flow analysis ||^). As a consequence, the 
optimal weight for a given particle is its flow Vn{pT,y) itself p3[ . Thus, a thorough flow analysis should go twice 



through Sees. |l| and [II. The first time, integrated flow can be extracted using some reasonable guess for the weights, 
thereby obtaining values for Vn{pT,y)', in turn, these values will serve as weights in the second, final analysis [ p^ . 

In Eq. (^), the number M of particles should be the same for all events: in each event, a set of M particles must 
be randomly chosen out of the Mtot detected particles. 

In order to obtain the cumulants, one first averages G'„(z) over events, which yields an average generating function 
{G,,{z)). We then define (w) 

C„(z) = m[(G„(z))^/*'-i] . (4) 

The cumulant of 2fc-particle correlations c„{2A;} is the coeflicient of z^z*^ j (/c!)^ in the power-series expansion of C„(2:). 
To construct the first three cumulants, one may truncate the series to order |z|^ and compute C„(z) at the following 
interpolation points: 

Zp,q = Xp^q + iyp^q, Xp^q = r^y/p COS ( j/p,, = tq sin ( j , (5) 

\ ^max / \ ^max / 

for p = 1,2,3 and (7 = 0,..., gmax ~ 1, where (/max > 8. The parameter tq must be chosen as a compromise between 
errors due to higher order terms in the power-series expansion, which rapidly increase with tq, and numerical errors. 
Assuming that the numerical error is proportional to the total number of elementary operations performed (of order 

MiVevts), we obtain the estimate ro ~ (tN^l'^^MY^^ M / (w^) where e is the accuracy of elementary operations, 
typically 10~^^ in double precision. This gives rp ~ 2 with weights of order unity and standard values of the 
multiplicity {M ^ 300) and number of events (iVcvts ~ 20000). 

Once the values Cn{zp.q) have been computed, they must be averaged over the phase of z: 

Cp= P= 1,2,3. (6) 

gmax 
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The cumulants of 2-, 4- and 6-particle correlations are then given respectively by 



c„{2}-- 3C1--C2 + -C3 , c„{4} = -(-5C1+4C2- C3), c„{6} = ^(3Ci-3C2 + C3). (7) 



These cumulants are related to the weighted integrated flow V„ = (we™'^'^"*^^). This is the point where acceptance 
considerations come into play. If the detector has full azimuthal coverage, each cumulant c„{2fc} gives an estimate of 
the corresponding Vn, which we denote by Vn{2k}: 

F„{2}2 = c„{2}, K{4}^ = -c„{4}, l/„{6}6 = c„{6}/4. (8) 



Generalized relations valid for detectors with partial azimuthal coverage or efficiency are given in Sec. [V . 

In practice, the use of higher order cumulants is often limited by statistics. The order of magnitude of statistical 
errors can easily be estimated. The computation of the cumulant c„{2fc} relies on the choice of 2k particles among M 
in each of the A'evts available events, i.e., it involves roughly AP^N^vts (2fc)-uplets of particles. Taking into account 

the weights w, the resulting statistical uncertainty on c„{2fc} is of order ©((w^)'^ / i/M^'^TVevts)- The relative error 
on Vn{2k} is thus of order 



6Vn{2k} 1_ 

Vn{2k} ^ xl'VK- 



(9) 



with Xn = Vn\/M/ {w'^). A thorough calculation, performed in Ref. @, Appendix D, shows that this order of magni- 
tude is indeed correct as long as x„ is not larger than unity, which is the case for most experiments at ultrarelativistic 
energies. 

In order to increase the statistics, it is possible to combine different runs performed in a given experiment, each of 
which has its own characteristics (different orientation of the magnetic field, etc.). For each of these runs, following the 
procedure described above leads to cumulants Cn^a{2k} (where a labels the run) and finally, accounting for the specific 
acceptance corrections, to flow estimates Vn,a{2k}^'^ . The proper way to combine the runs consists in averaging the 
Vn,a{2k}^'' {not the Vn,a{2k}) of the various runs, weighted by the number of events in each run. Note that if 
statistical fluctuations are large, it may happen that a run yield a negative value of Vn^a{2k}^'' . Such a run must 
nevertheless be included in the averaging procedure. 



III. DIFFERENTIAL FLOW 



Let us turn to differential flow, i.e., to the flow of an identified particle in a restricted portion of phase space 
(typically a narrow px or y interval). We call such a particle a "proton", and denote its azimuth by ip] thus, the 
differential fiow is v' = (e™(^-*«)), where the average value is taken over all protons. 

Experimentally, the differential flow is obtained from the azimuthal correlations between the proton and the particles 
previously used to determine the integrated flow, which we call "pions" . It is well known in the standard flow analysis 
that from the reaction plane reconstructed in harmonic n, one may reconstruct not only the corresponding v'^, but 
also higher harmonics w^j„, where m is an integer (in practice, m = 1 or 2). Here, this is done by correlating the 
proton with m pions, i.e., by measuring (exp[in(mV' — (f>i — ■ ■ ■ — 0m)])- As in the case of integrated flow, nonflow 
contributions to this correlation can be eliminated by going to higher orders, i.e., by correlating the proton with 
2k + m pions with fc > (in practice fc = or 1), and constructing a cumulant dmn/n{2k + m + 1}. The subscript 
refers to the fact that is measured by using pions in harmonic n, while the number in curly brackets is the order 
of the correlation: 2fc + m pions and 1 proton. 

This cumulant retains only the contributions from flow, proportional to v'^^V^^^"^ , and from {2k-\-m-\- l)-particle 
nonflow correlations, which scale like 0(1/Af^*^+'"): all lower order nonflow correlations have been removed. If the 
contribution of flow dominates, the cumulant dmn/n{2k + m + 1} yields an estimate of which we naturally denote 

V' , |2fc + TO + 1|. 

As in the case of integrated flow, we derive the cumulants (im„/„{2fc + m + 1} using a generating function, namely 



(e™"'^ G„(z)) 
(Gniz)) 



T^mn/niz) = ■ (10) 



The cumulant dmn/n{2k+m+l} is the real part of the coefficient of z*'' z'^^"'' /[k\{k+m)\] in the power-series expansion 
of 2?^„/„(z). In the numerator of Eq. (|lO|), the average is performed over all protons (i.e., an event with 2 protons is 
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counted twice; this was not stated correctly in On the other hand, the denominator is averaged over all events. 

If the proton is one of the "pions", i.e., if it was used in the calculation of the generating function Eq. (^), one should 
divide Gn{z) by 1 + Wj{z*e"'''^ + ze~™'^)/M, where -0 is the proton azimuth, to avoid autocorrelations. Note that 
while the number of pions in Eq. was fixed, the number of protons must be allowed to fluctuate from event to 
event: to increase statistics, one should use all available protons. 

To extract the cumulants, one computes the product 2;*'"I?m„/„(z) at the points Zp^q, Eq. (||); one then takes the 
real part, and averages over angles: 



cos I m I Xp^q + sin I m I Yp q 



(11) 



with p — 1,2,3 and Xp^q + iYp ^ = I'mn/n(^p.g)- Note that although we present the integrated and differential flow 
analyses as two successive steps, they are in fact simultaneous: while the generating function Gn(z) is calculated for 
a given event, one can at the same time compute its product by e™"''' for the numerator of Eq. (|l(]). 
For TO = 1 (useful for v[^^ or ^2/2)1 the lowest order cumulants are given by 

d„/„{2} = l(2Di-ii^2), d„/„{4} = ^(-2Di + i?2), (12) 
while for to = 2, which can be used to derive Wj/i' 

d2„/„{3} = ;i (4 ^1 - ^ ' d^n/ni^} = ^ ^1 + ^ . (13) 

These cumulants must then be related to the differential flow. For a perfect detector: 

</„{2} = d„/«{2}/K, </„{4} = -d,,/„{4}/V^, 

V2n/n{^} = d2n/nm/Vl v',^/„ {5} ^ - d,,,/,,{5} / i2V^) . (14) 

Note that these relations involve the integrated flow Vn obtained in the previous section. Relations for a nonisotropic 
acceptance are given in next Section. The order of magnitude of the statistical uncertainty can be estimated using 
the same arguments as for integrated flow. One obtains 

SvLn/n{2k + m + l}^ , (15) 

Xn V J V 

where N' is the number of protons used in the analysis. It may be worth noticing that the statistical error does not 
only depends on the available event statistics and multiplicity. It also strongly depends on the flow itself through the 
parameter Xn- 



IV. ACCEPTANCE CORRECTIONS 



Our method can be used even when the detector used to measure the particle azimuths has only partial azimuthal 
coverage, provided that the event centralities be determined with an independent detector, as, e.g., a ZDC, with an 
approximately isotropic coverage. This is to make sure that the apparent multiplicity/centrality is not biased by the 
orientation of the reaction plane with respect to the detector. 

A specific detector is characterized by its acceptance/efficiency function A{j,(j>,pT,y), which represents the prob- 
ability that a particle of type j (pion, proton, etc.) with azimuth 0, transverse momentum px, and rapidity y, be 
detected. In practice, A{j, 0,pt, y) is proportional to the number of hits in a {4>,pt, y) bin, and thus can be obtained 
in a straightforward way while scanning through the data. The Fourier coefficients of the acceptance function are 

A,,{j,pT,y)= f A{j,q^,pT,y)e-'^U<j). (16a) 
Jo 

These differential coefficients can be integrated, with appropriate weighting and a sum over the various types of 
particles used for the flow analysis, so as to describe the "integrated" acceptance of the detector: 
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j 

/ w{j,pT,y)Ao{j,pT,y)dpTdy 



(16b) 



Note our introducing the weights w{j,pT,y), which are of course the same as in Eq. (|^). 

When the acceptance of a detector is not perfectly isotropic in (j), the cumulants Eq. (|^) wiU mix different flow 
harmonics: Cn{2k} does not depend only on Vn, but also on other Vp with p ^ n. For instance, 



ci{4} = - [(1 - |ai|')* + 4(1 - |ai|')'|a2 - ajf + \a2 - a?|*] V* 

- [|ai - a2al\'* + 4 |ai - a2al\'^\a3 - aia2p + jas - aia2|'*] V2 



and 



C2{4} 



[lal 

[(1- 



' a2ai| 



' 4 \ai 
4(1- 



a2a^|'|a3 - aia2|' + [03 - aia2|'*] V-^ 



|a2n'|a4 ■ 



,2|2 



a4 



,2|41 
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2 J 



(17a) 



(17b) 



where we have assumed that all other flow harmonics Vn>3 are negligible. The corresponding relations for the 
cumulants ci{2} and C2{2} can be found in Ref. jl4|, Appendix CI. If the detector acceptance is not too bad, the 
coefficients afc^p will be small, and these expressions are close to Eqs. (||), valid for perfect detector. 

Equations ( pT[ ) form a linear system which can easily be inverted to express Vi and V2 (or, more precisely, the 
estimates Vi{4}^ and V2{4}^) as functions of ci{4} and C2{4}. 

Let us now consider differential flow. Integrated and differential flows may be measured using two different de- 
tectors: for instance, a large acceptance detector for integrated flow, and a smaller one, but with better particle 
identification or px determination, for differential flow. For sake of generality, we thus denot e by A'{j,il),px,y) the 
corresponding acceptance function and by A'^{j,pT,y) its Fourier coefficients defined as in (16a). The differential 
acceptance coefficients a'j, are then defined as in Eq. ( |16b| ) , without the weights and the summation over j (since one 
usually measures the differential flow of identified particles) and with the integration over px and y restricted to the 
phase-space region under interest (typically, one integrates over pT or y, so as to obtain Vn as a function of y or px, 
respectively) . 

Once again, anisotropics in the detector acceptance lead to some interference between the fiow harmonics in the 
expressions of the cumulants dmn/n{'^k + m + 1}. Here are, for example, the relations between the lowest order 
cumulants and flow which allow one to extract the differential directed flow v'l and the differential elliptic flow V2 
obtained either with respect to Vi or V2, i.e., what we denote v^-^ and (relations for other cumulants can be 
found in Ref. |l|. Appendix C2): 



rfi/i{4} 



-Re 
-Re 



rf2/l{3} = 

rf2/2{4} - -Re 
-Re 



(1 - |ai|2) ((1 - |ai|2)2 + 2 |a2 - a?|') + (4)*(a2 - a?) (2(1 - \ai\Y + \o2 - a?| 
a[{al — a^ai) ^|a* — a^oif -|- 2 jas — aia2|'^ 
+ {a'^)*{a3 - aia2) (2 \al - a2ai|' -|- jaa - aia2|'^ i'2^2^: 
Re[{l-\a,ff + {aina2-alf]v'2V,^ 

-I- 2 Re [ia[)*{al - alai){a2 - al) + {a^yias - a2ai)(l - |ai|^)] v[ViV2, 

(1 - laaH ((1 - laal')' + 2 [04 - a^f ) + (a'^yia^ - a^) (2(1 - jasH^ + - aj] 
{a[y{ai - a2al) (^|a* - ajai]' -|- 2 |a3 - aia2|'^ 
+ {a'^y{a3 - aia2) (2 \al - ajai]' + |a3 - aia2n v[Vi . 



V'2V^ 



(18a) 
(18b) 

(18c) 



In these expressions, we have neglected terms involving higher flow harmonics Vn>3 or f^>3, and Re means that one 
should take the real part. Of course, they reduce to Eqs. (^J) when the acceptance is isotropic. 

Any two of Eqs. (^ constitute a linear system which can be inverted to obtain v'-^ and V2 once Vi and V2 have 
been extracted. 
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V. ADDITIONAL REMARKS AND COMMENTS 



The cumulant method not only eUminates nonflow eorrelations, it also provides several independent estimates of 
the flow from cumulants of various orders. Since nonflow correlations between 4 or more particles are expected to 
be negligible, Ki{4} and Ki{6} should be consistent with each other within statistical error bars. With the high 
multiplicity produced at ultrarelativistic energies and the large acceptance detectors available, estimates show that 
one could even construct cumulants of 8-, 10-particle correlations. Since the generating function formalism yields all 
cumulant orders at once, one would simply need to increase the number of interpolation points Zp^q, which would 
result in a moderate increase in computer time. Checking that all higher order cumulants yield compatible values of 
Vn would give the first direct evidence that azimuthal correlations are of collective origin. Similarly, in the case of 
differential flow, various cumulants orders yield independent estimates which can be compared to one another. 

Other consistency checks can be proposed to test the reliability of the results. First, one should perform the analysis 
with at least two values of the parameter rg entering Eq. (^), in order to check the stability against numerical errors. 
Then, one can try different values of M in Eq. (||), as, for instance, M = 0.8 (Aftot) and M = 0.6 (A/tot)- That may 
also be a way to increase the statistics. In the same spirit, it is still possible to let M vary from event to event in 
Eq. (^), taking M — Mt^t, especially if the acceptance is reasonably good. In that case, M must be replaced by the 
average (A/tot) in Eq- (H), and one should check that the results are consistent with what would be obtained with a 
fixed M. 

The last consistency check regards differential flow. We have seen that the integrated flow Vn can be calculated 
using different weights wj; the resulting weighted averages Vn's differ. However, the values of the differential flow 
t)^^^^{2fc + m+ 1} obtained with different weights should be consistent within error bars if they are not contaminated 
by nonflow effects. 

The recent STAR analysis |l5) illustrates quite well the relevance of cumulants to the flow analysis: for integrated 
flow, the lowest order y„{2} reproduces the results of the standard two-particle methods, as it should. On the other 
hand, the higher order estimate Ki{4} is lower than y„{2} beyond statistical error bars, especially for peripheral 
collisions, thereby suggesting that nonflow correlations are important. Although statistical uncertainties on higher 
order cumulants are larger, this loss is therefore compensated by the gain on errors from nonflow effects. 
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